Setup

library(seriation)
library(dplyr)
library(purrr)
library(ggplot2)
library(phyloseq)
library(microViz)
library(shiny)
mice <- readRDS(here::here('data/mice.rds'))
load(here::here('data/shao19.rda'))
shao19 <- tax_fix(shao19, verbose = FALSE) %>% tax_names2rank("species")
shao4 <- shao19 %>% ps_filter(family_role == 'child', infant_age == 4) 

Beta diversity (part 2)

Euclidean distances

What about Euclidean distances? What are those?

  • Show pythagoras \(c = \sqrt{a^2 + b^2}\) and explain generalisation to more dimensions, where each dimension is the counts of one taxon.

Issues

  • Sensitive to sparsity (double-zero problem) –> filter rare taxa
  • Excessive emphasis on high-abundance taxa –> transform features first
  • The PCoA looks weird! most samples bunched in the middle with spindly projections..
shao4 %>% 
  tax_filter(min_prevalence = 2.5/100, verbose = FALSE) %>% 
  tax_agg(rank = 'genus') %>% 
  dist_calc(dist = 'euclidean') %>% 
  ord_calc(method = 'PCoA') %>% 
  ord_plot(alpha = 0.6, size = 2) +
  theme_classic(12) +
  coord_fixed(0.7) + 
  geom_rug(alpha = 0.1)

Abundance transformation

We already did two transformations with tax_transform(): binary (for Binary Jaccard distances) and compositional (for barplots).

Now we need log transformations, and the centered-log-ratio, CLR, transformation.

Log transformation

First let’s look at the abundance again, this time with heatmaps.

# Getting the taxa in abundance order up front
# to keep it consistent across multiple plots
shao4_sorted <- shao4 %>% 
  tax_sort(by = sum, at = 'genus', trans = 'compositional', tree_warn = FALSE)

Each column is a sample (from an infant), and each row is a taxon.

shao4_sorted %>% 
  tax_transform(trans = 'identity', rank = 'genus') %>% 
  comp_heatmap(
    samples = 1:30, taxa = 1:20, grid_lwd = 2, name = 'Counts',
    tax_seriation = 'Identity', sample_seriation = "Identity"
  )

shao4_sorted %>% 
  tax_transform(trans = 'compositional', rank = 'genus') %>% 
  comp_heatmap(
    samples = 1:30, taxa = 1:20, grid_lwd = 2, name = 'Prop.',
    tax_seriation = 'Identity', sample_seriation = "Identity"
  )

Even though we have picked the top 20 most abundant genera, there are still a lot of zeros, We need to deal with the zeros, because log(0) is undefined. The solution is to add a small amount to every value (or just every zero), before applying the log transformation. This small value is often called a pseudo-count.

What value should we use for the pseudo-count?

One option is to just add 1, and another popular option is to add half of the smallest observed real value (from across the whole dataset).

shao4_sorted %>% 
  tax_transform(rank= 'genus', trans = 'log10', zero_replace = 1) %>% 
  comp_heatmap(
    samples = 1:30, taxa = 1:20, grid_lwd = 2, name = 'log10\n(x+1)', 
    tax_seriation = 'Identity', sample_seriation = "Identity"
  )

shao4_sorted %>% 
  tax_agg(rank = 'genus') %>% 
  # tax_transform(trans = 'compositional') %>% # compositional also possible
  tax_transform(trans = 'log10', zero_replace = 'halfmin', chain = TRUE) %>% 
  comp_heatmap(
    samples = 1:30, taxa = 1:20, grid_lwd = 2, name = 'log10\nhalfmin', 
    tax_seriation = 'Identity', sample_seriation = "Identity"
  )

Log-transformation, zero replacement summary (keep it simple and record your approach).

Centered log ratio transformation:

Compositionality - centered-log-ratio transformation

The centered log-ratio (clr) transformation uses the geometric mean of the sample vector as the reference

shao4_sorted %>% 
  tax_agg(rank = 'genus') %>% 
  # tax_transform(trans = 'compositional') %>% # compositional also possible
  tax_transform(trans = 'clr', zero_replace = 'halfmin', chain = TRUE) %>% 
  comp_heatmap(
    samples = 1:30, taxa = 1:20, grid_lwd = 2, name = 'CLR\nhalfmin', 
    colors = heat_palette(sym = TRUE),
    tax_seriation = 'Identity', sample_seriation = "Identity"
  )

Overview of CoDa problem - slides only.

The sequencing data gives us relative abundances, not absolute abundances. The total number of reads sequenced per sample is an arbitrary total.

If one taxon blooms, whilst everything else stays stable, the relative abundance of all other taxa must (appear to) go down.

This leads to two main types of problem:

  • interpretation caveats: see differential abundance section later
  • statistical issues: taxon abundances are not independent, but (weakly?) negatively correlated

This is worse with simpler ecosystems. There is the same problem in theory with RNAseq data, but I suspect it is less bothersome because there are many more competing “species” of RNA transcript than there are bacterial species in even a very complex microbiome.

The centered-log-ratio transformation (along with some other similar ratio transformations) are claimed to help with the statistical issues by transforming the abundances from the simplex to the real space.

Practically, the CLR transformation involves finding the geometric mean of each sample, and then dividing abundance of each taxon in that sample by this geometric mean. Finally you take the natural log of this ratio.

For more details, check out Gloor 2017, and work by Thomas Quinn. TODO: link.

PCA

Principal Components Analysis.

Quite similar to Principal Co-ordinates Analysis.

In fact, PCA produces equivalent results to PCoA with euclidean distances. So let’s perform the CLR-transform first and check PCA and euclidean PCoA are the same.

shao4 %>% 
  tax_filter(min_prevalence = 2.5/100, verbose = FALSE) %>% 
  tax_transform(rank = 'genus', trans = 'clr', zero_replace = 'halfmin') %>% 
  dist_calc(dist = 'euclidean') %>% 
  ord_calc(method = 'PCoA') %>% 
  ord_plot(alpha = 0.6, size = 2, color = 'birth_mode') +
  theme_classic(12) +
  coord_fixed(0.7) 

shao4 %>% 
  tax_filter(min_prevalence = 2.5/100, verbose = FALSE) %>% 
  tax_transform(rank = 'genus', trans = 'clr', zero_replace = 'halfmin') %>% 
  ord_calc(method = 'PCA') %>% 
  ord_plot(alpha = 0.6, size = 2, color = 'birth_mode') +
  theme_classic(12) +
  coord_fixed(0.7)

So why is PCA interesting for us? Because the Principal components are built directly from a (linear) combination of the original features.

That means we know how much each taxon contributes to each PC axis, and we can plot this information (loadings) as arrows, alongside the sample points.

pca <- shao4 %>% 
  tax_filter(min_prevalence = 2.5/100, verbose = FALSE) %>% 
  tax_transform(rank = 'genus', trans = 'clr', zero_replace = 'halfmin') %>% 
  ord_calc(method = 'PCA') %>% 
  ord_plot(
    alpha = 0.6, size = 2, color = 'birth_mode',
    plot_taxa = 1:6, tax_vec_length = 0.275,
    tax_lab_style = tax_lab_style(
      type = 'text', max_angle = 90, aspect_ratio = 0.7, 
      size = 3, fontface = 'bold'
    ),
  ) +
  theme_classic(12) +
  coord_fixed(0.7, clip = 'off')
pca

How to interpret the taxa loading vectors? Cautiously.

The relative length and direction of an arrow indicates how much that taxon contributes to the variation on each visible PC axis. e.g. Variation in Enterococcus contributes quite a lot to variation along the PC2 axis.

This allows you to infer that samples positioned at the bottom of the plot will tend to have higher relative abundance of Enterococcus than samples at the top of the plot.

Interestingly, samples on the right of the plot (which tend to be vaginally-delivered infants) seem to have relatively more Bifidobacterium, Bacteroides and Escherichia, whilst the C-section born infants have relatively more Veillonella.

There are caveats and nuance to the interpretation of these plots, which are called PCA bi-plots, and you can read more about here: TODO link to appropriate gusta me ecology page.

(Side note, Phocaeicola were considered part of Bacteroides until this year!)

Iris plot

You might have already noticed this pattern, when exploring and making barplots interactively with ord_explore earlier. We can make another kind of barplot now, using the PCA information to order our samples in a circular layout.

iris <- shao4 %>% 
  tax_filter(min_prevalence = 2.5/100, verbose = FALSE) %>% 
  tax_transform(rank = 'genus', trans = 'clr', zero_replace = 'halfmin') %>% 
  ord_calc(method = 'PCA') %>% 
  ord_plot_iris(
    tax_level = 'genus', n_taxa = 12, other = "Other", 
    anno_colour = 'birth_mode', anno_colour_style = list(alpha = 0.6, size = 0.6, show.legend = FALSE)
  )
iris

patchwork::wrap_plots(pca, iris, nrow = 1, guides = 'collect')

Taxon stats

From the PCA loadings and barplots above, we have some strong suspicions about which taxa have a higher relative abundance in vaginally delivered infants than in c-section delivered infants, and vice versa, but we can also statistically test this. This is often called “differential abundance” testing, in the style of “differential expression” testing from the transcriptomics field.

shao4 %>% 
  comp_barplot(
    tax_level = 'genus', n_taxa = 12, facet_by = 'birth_mode',
    label = NULL, bar_outline_colour = NA
  ) +
  coord_flip() +
  theme(axis.ticks.y = element_blank())

Model one taxon

We will start by creating a linear regression model for one genus, Bacteroides. We will transform the count data by first making it proportions, and then taking the binary logarithm, log2, after adding a pseudocount.

bacteroidesRegression1 <- shao4 %>% 
  tax_transform("compositional", rank = "genus") %>% 
  tax_transform("log2", zero_replace = "halfmin", chain = TRUE) %>%
  tax_model(type = "lm", rank = "genus", taxa = "Bacteroides", variables = "birth_mode") %>% 
  pluck(1) 
## Modelling: Bacteroides
# looking at the regression results
summary(bacteroidesRegression1)
## 
## Call:
## Bacteroides ~ birth_mode
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.7492 -0.6172 -0.6172  2.6421 18.0804 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -19.3756     0.4863  -39.84   <2e-16 ***
## birth_modevaginal   7.1320     0.6812   10.47   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.957 on 304 degrees of freedom
## Multiple R-squared:  0.265,  Adjusted R-squared:  0.2626 
## F-statistic: 109.6 on 1 and 304 DF,  p-value: < 2.2e-16
confint(bacteroidesRegression1)
##                        2.5 %     97.5 %
## (Intercept)       -20.332614 -18.418542
## birth_modevaginal   5.791662   8.472414
broom::tidy(bacteroidesRegression1, conf.int = TRUE)
## # A tibble: 2 × 7
##   term              estimate std.error statistic   p.value conf.low conf.high
##   <chr>                <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
## 1 (Intercept)         -19.4      0.486     -39.8 1.08e-122   -20.3     -18.4 
## 2 birth_modevaginal     7.13     0.681      10.5 4.13e- 22     5.79      8.47

Click here for optional ggplot2 extension exercise:

Starting from a dataframe like the one produced by the code below, plot:

  1. Easy: The percentage prevalence of Bacteroides in each birth_mode group
  2. Medium: The distribution of relative abundance of Bacteroides in each birth_mode group, omitting zeros, on a log2 scale
  3. Hard: Do task 1 or 2 for for several taxa in one plot - (hint: pivot_longer)
  shao4 %>% 
  tax_transform("compositional", rank = "genus") %>% 
  ps_get() %>% 
  ps_otu2samdat(taxa = 'Bacteroides') %>% 
  samdat_tbl()

We can fit a model with covariates, as we did for PERMANOVA. We are going to convert the categorical variables into indicator (dummy) variables, and scale the continuous covariates to 0 mean and SD 1 (z-scores). You’ll see this will make our subsequent plots easier to interpret later.

shao4 <- shao4 %>% 
  ps_mutate(
    C_section = if_else(birth_mode == 'c_section', true = 1, false = 0),
    Female = if_else(sex == 'female', true = 1, false = 0),
    Birth_weight_Z = scale(birth_weight, center = TRUE, scale = TRUE),
    Reads_Z = scale(number_reads, center = TRUE, scale = TRUE)
  )
bacteroidesRegression2 <- shao4 %>% 
  tax_transform("compositional", rank = "genus") %>% 
  tax_transform("log2", zero_replace = "halfmin", chain = TRUE) %>%
  tax_model(
    type = "lm", rank = "genus", taxa = "Bacteroides", 
    variables = c('C_section', 'Female', 'Birth_weight_Z', 'Reads_Z')
  ) %>% 
  pluck(1) 
## Modelling: Bacteroides
# looking at the regression results
summary(bacteroidesRegression2)
## 
## Call:
## Bacteroides ~ C_section + Female + Birth_weight_Z + Reads_Z
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.4271 -2.1555 -0.4115  2.8176 18.1784 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -11.7942     0.6103 -19.325   <2e-16 ***
## C_section       -7.5696     0.7206 -10.505   <2e-16 ***
## Female          -0.3809     0.7101  -0.536    0.592    
## Birth_weight_Z   0.3277     0.3514   0.932    0.352    
## Reads_Z          0.5361     0.3620   1.481    0.140    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.934 on 286 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.2854, Adjusted R-squared:  0.2754 
## F-statistic: 28.55 on 4 and 286 DF,  p-value: < 2.2e-16
broom::tidy(bacteroidesRegression2, conf.int = TRUE)
## # A tibble: 5 × 7
##   term           estimate std.error statistic  p.value conf.low conf.high
##   <chr>             <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 (Intercept)     -11.8       0.610   -19.3   8.15e-54  -13.0      -10.6 
## 2 C_section        -7.57      0.721   -10.5   4.81e-22   -8.99      -6.15
## 3 Female           -0.381     0.710    -0.536 5.92e- 1   -1.78       1.02
## 4 Birth_weight_Z    0.328     0.351     0.932 3.52e- 1   -0.364      1.02
## 5 Reads_Z           0.536     0.362     1.48  1.40e- 1   -0.176      1.25

Many methods

This method is what MaAsLin2 does by default (except they call the compositional transformation “Total Sum Scaling (TSS)”). This is quite an uncomplicated method, so we will stick with this for today, but many statistical methods have been developed for differential abundance analyses.

Microbiome abundance data are quite awkward, statistically speaking, due to their sparseness and compositionality. Each successive method claims to handle some aspect of this awkwardness “better” than previous methods.

The aim is to have a method with adequate power to detect true associations, whilst controlling the type 1 error rate, the “false positive” reporting of associations that are not “truly” present.

Results are surprisingly inconsistent across the different methods, as demonstrated this year in a fascinating analysis by Jacob Nearing and colleagues.

What to do?

  • Filter out the noise & interpret results with caution! use multiple testing corrections
  • Remember it’s all relative (abundance)
  • Try multiple methods and/or use same method as previous study if replicating - Avoid Lefse and edgeR? - Beware: Not all methods allow covariate adjustment & few allow random effects (for time-series)

Model all the taxa!

We’re not normally interested in just one taxon! And often it’s also hard to decide which taxonomic rank we are most interested in modelling!

Lower ranks like species or ASVs give better resolution but also more sparsity and classification uncertainty… Higher ranks e.g. classes, could also be more powerful if you think most taxa within that class will follow a similar pattern.

Insert meme here?

So now we will fit a similar model for almost* every taxon at every rank we have available, from phylum down to species.

*We’ll actually filter out species with a prevalence of less than 10%.

# The code for `taxatree_models` is quite similar to tax_model. 
# However, you might need to run `tax_prepend_ranks` to ensure that each taxon at each rank is always unique. 
shaoModels <- shao4 %>% 
  tax_prepend_ranks() %>%
  tax_transform("compositional", rank = "species", keep_counts = TRUE) %>% 
  tax_filter(min_prevalence = 0.1, undetected = 0, use_counts = TRUE) %>% 
  tax_transform(trans = "log2", chain = TRUE, zero_replace = "halfmin") %>% 
  taxatree_models(
    type = lm, 
    ranks = c("phylum", "class", "order", "family", "genus", "species"), 
    variables = c('C_section', 'Female', 'Birth_weight_Z', 'Reads_Z')
  )
## Proportional min_prevalence given: 0.1 --> min 31/306 samples.
## 2022-05-27 17:55:56 - modelling at rank: phylum
## 2022-05-27 17:55:56 - modelling at rank: class
## 2022-05-27 17:55:56 - modelling at rank: order
## 2022-05-27 17:55:56 - modelling at rank: family
## 2022-05-27 17:55:56 - modelling at rank: genus
## 2022-05-27 17:55:57 - modelling at rank: species
shaoModels
## ps_extra object - a list with phyloseq and extras:
## 
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 39 taxa and 306 samples ]
## sample_data() Sample Data:       [ 306 samples by 15 sample variables ]
## tax_table()   Taxonomy Table:    [ 39 taxa by 6 taxonomic ranks ]
## 
## ps_extra info:
## tax_agg = species tax_transform = compositional&log2
## 
## $counts OTU Table: [ 39 taxa and 306 samples ]
## 
## $taxatree_models list:
## Ranks: phylum/class/order/family/genus/species

Why filter the taxa? It’s less likely that we are interested in rare taxa, and models of rare taxon abundances are more likely to be unreliable. Reducing the the number of taxa modelled also makes the process faster and makes visualizing the results easier!

Getting stats from the models

Next we will get a data.frame containing the regression coefficient estimates, test statistics and corresponding p values from all these regression models.

shaoStats <- taxatree_models2stats(shaoModels)
shaoStats
## ps_extra object - a list with phyloseq and extras:
## 
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 39 taxa and 306 samples ]
## sample_data() Sample Data:       [ 306 samples by 15 sample variables ]
## tax_table()   Taxonomy Table:    [ 39 taxa by 6 taxonomic ranks ]
## 
## ps_extra info:
## tax_agg = species tax_transform = compositional&log2
## 
## $counts OTU Table: [ 39 taxa and 306 samples ]
## 
## $taxatree_stats dataframe:
## 86 taxa at 6 ranks: phylum, class, order, family, genus, species 
## 4 terms: C_section, Female, Birth_weight_Z, Reads_Z
shaoStats$taxatree_stats
## # A tibble: 344 × 7
##    term           taxon             rank   estimate std.error statistic  p.value
##    <fct>          <chr>             <fct>     <dbl>     <dbl>     <dbl>    <dbl>
##  1 C_section      p: Proteobacteria phylum  -3.44       1.31    -2.62   9.23e- 3
##  2 Female         p: Proteobacteria phylum   0.595      1.29     0.460  6.46e- 1
##  3 Birth_weight_Z p: Proteobacteria phylum   0.0179     0.640    0.0279 9.78e- 1
##  4 Reads_Z        p: Proteobacteria phylum  -1.22       0.659   -1.86   6.44e- 2
##  5 C_section      p: Actinobacteria phylum -16.9        2.37    -7.12   8.98e-12
##  6 Female         p: Actinobacteria phylum  -2.80       2.33    -1.20   2.30e- 1
##  7 Birth_weight_Z p: Actinobacteria phylum   1.65       1.15     1.43   1.55e- 1
##  8 Reads_Z        p: Actinobacteria phylum  -2.86       1.19    -2.41   1.67e- 2
##  9 C_section      p: Firmicutes     phylum  15.3        4.21     3.62   3.43e- 4
## 10 Female         p: Firmicutes     phylum   0.741      4.15     0.179  8.58e- 1
## # … with 334 more rows

Adjusting p values

As we have performed a lot of statistical tests here, it is quite possible that could we find some significant p-values by chance alone.

So we should correct for multiple testing / control the false discovery rate or family-wise error rate.

Instead of applying these adjustment methods across all 86 taxa models at all ranks, the default behaviour is to control the family-wise error rate per taxonomic rank.

shaoStats <- shaoStats %>% taxatree_stats_p_adjust(method = "BH", grouping = "rank")
# notice the new variable
shaoStats$taxatree_stats
## # A tibble: 344 × 8
## # Groups:   rank [6]
##    term          taxon rank  estimate std.error statistic  p.value p.adj.BH.rank
##    <fct>         <chr> <fct>    <dbl>     <dbl>     <dbl>    <dbl>         <dbl>
##  1 C_section     p: P… phyl…  -3.44       1.31    -2.62   9.23e- 3      2.95e- 2
##  2 Female        p: P… phyl…   0.595      1.29     0.460  6.46e- 1      7.38e- 1
##  3 Birth_weight… p: P… phyl…   0.0179     0.640    0.0279 9.78e- 1      9.78e- 1
##  4 Reads_Z       p: P… phyl…  -1.22       0.659   -1.86   6.44e- 2      1.47e- 1
##  5 C_section     p: A… phyl… -16.9        2.37    -7.12   8.98e-12      7.19e-11
##  6 Female        p: A… phyl…  -2.80       2.33    -1.20   2.30e- 1      3.07e- 1
##  7 Birth_weight… p: A… phyl…   1.65       1.15     1.43   1.55e- 1      2.48e- 1
##  8 Reads_Z       p: A… phyl…  -2.86       1.19    -2.41   1.67e- 2      4.46e- 2
##  9 C_section     p: F… phyl…  15.3        4.21     3.62   3.43e- 4      1.83e- 3
## 10 Female        p: F… phyl…   0.741      4.15     0.179  8.58e- 1      9.16e- 1
## # … with 334 more rows

Plot all the taxatree_stats!

taxatree_plots() allows you to plot statistics (e.g. point estimates and significance) from all of the taxa models onto a tree layout. The taxon models are organised by rank, radiating out from the central root node from e.g. Phyla around the center to Species in the outermost ring.

taxatree_plots() itself returns a list of plots, which you can arrange into one figure with the patchwork package for example (and/or cowplot).

shaoStats %>% 
  taxatree_plots(node_size_range = c(1, 3), sig_stat = 'p.adj.BH.rank') %>% 
  patchwork::wrap_plots(ncol = 2, guides = "collect")

Taxatree Key

But how do we know which taxa are which nodes? We can create a labelled grey tree with taxatree_plotkey(). This labels only some of the taxa based on certain conditions that we specify.

set.seed(123) # label position 
key <- shaoStats %>% 
  taxatree_plotkey(
    taxon_renamer = function(x) stringr::str_remove(x, "[pfg]: "),
    # conditions below, for filtering taxa to be labelled
    rank == "phylum" | rank == "genus" & prevalence > 0.2
    # all phyla are labelled, and all genera with a prevalence of over 0.2
  ) 
key

You can do more with these trees to customise them to your liking. See an extended tutorial here on the microViz website: including how to directly label taxa on the colored plots, change the layout and style of the trees, and even how to use a different regression modelling approach.

Fin.

End workshop with recap and notes/pointers on current trends and topics not covered?

Session info

Records your package versions etc. Useful for debugging / reproducing analysis.

devtools::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.2.0 (2022-04-22)
##  os       macOS Monterey 12.4
##  system   aarch64, darwin20
##  ui       X11
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       Europe/Amsterdam
##  date     2022-05-27
##  pandoc   2.17.1.1 @ /Applications/RStudio.app/Contents/MacOS/quarto/bin/ (via rmarkdown)
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  ! package          * version  date (UTC) lib source
##    ade4               1.7-19   2022-04-19 [2] CRAN (R 4.2.0)
##    ape                5.6-2    2022-03-02 [2] CRAN (R 4.2.0)
##    assertthat         0.2.1    2019-03-21 [2] CRAN (R 4.2.0)
##    backports          1.4.1    2021-12-13 [2] CRAN (R 4.2.0)
##    Biobase            2.56.0   2022-04-26 [2] Bioconductor
##    BiocGenerics       0.42.0   2022-04-26 [2] Bioconductor
##    biomformat         1.24.0   2022-04-26 [2] Bioconductor
##    Biostrings         2.64.0   2022-04-26 [2] Bioconductor
##    bitops             1.0-7    2021-04-24 [2] CRAN (R 4.2.0)
##    brio               1.1.3    2021-11-30 [2] CRAN (R 4.2.0)
##    broom              0.8.0    2022-04-13 [2] CRAN (R 4.2.0)
##    bslib              0.3.1    2021-10-06 [2] CRAN (R 4.2.0)
##    cachem             1.0.6    2021-08-19 [2] CRAN (R 4.2.0)
##    Cairo              1.5-15   2022-03-16 [2] CRAN (R 4.2.0)
##    callr              3.7.0    2021-04-20 [2] CRAN (R 4.2.0)
##    circlize           0.4.15   2022-05-10 [2] CRAN (R 4.2.0)
##    cli                3.3.0    2022-04-25 [2] CRAN (R 4.2.0)
##    clue               0.3-60   2021-10-11 [2] CRAN (R 4.2.0)
##    cluster            2.1.3    2022-03-28 [2] CRAN (R 4.2.0)
##    codetools          0.2-18   2020-11-04 [2] CRAN (R 4.2.0)
##    colorspace         2.0-3    2022-02-21 [2] CRAN (R 4.2.0)
##    ComplexHeatmap     2.12.0   2022-04-26 [2] Bioconductor
##    corncob            0.2.0    2021-03-11 [2] CRAN (R 4.2.0)
##    crayon             1.5.1    2022-03-26 [2] CRAN (R 4.2.0)
##    data.table         1.14.2   2021-09-27 [2] CRAN (R 4.2.0)
##    DBI                1.1.2    2021-12-20 [2] CRAN (R 4.2.0)
##    desc               1.4.1    2022-03-06 [2] CRAN (R 4.2.0)
##    devtools           2.4.3    2021-11-30 [2] CRAN (R 4.2.0)
##    digest             0.6.29   2021-12-01 [2] CRAN (R 4.2.0)
##    doParallel         1.0.17   2022-02-07 [2] CRAN (R 4.2.0)
##    dplyr            * 1.0.9    2022-04-28 [2] CRAN (R 4.2.0)
##    ellipsis           0.3.2    2021-04-29 [2] CRAN (R 4.2.0)
##    evaluate           0.15     2022-02-18 [2] CRAN (R 4.2.0)
##    fansi              1.0.3    2022-03-24 [2] CRAN (R 4.2.0)
##    farver             2.1.0    2021-02-28 [2] CRAN (R 4.2.0)
##    fastmap            1.1.0    2021-01-25 [2] CRAN (R 4.2.0)
##    foreach            1.5.2    2022-02-02 [2] CRAN (R 4.2.0)
##    fs                 1.5.2    2021-12-08 [2] CRAN (R 4.2.0)
##    future             1.25.0   2022-04-24 [2] CRAN (R 4.2.0)
##    future.apply       1.9.0    2022-04-25 [2] CRAN (R 4.2.0)
##    generics           0.1.2    2022-01-31 [2] CRAN (R 4.2.0)
##    GenomeInfoDb       1.32.1   2022-04-28 [2] Bioconductor
##    GenomeInfoDbData   1.2.8    2022-05-11 [2] Bioconductor
##    GetoptLong         1.0.5    2020-12-15 [2] CRAN (R 4.2.0)
##    ggforce            0.3.3    2021-03-05 [2] CRAN (R 4.2.0)
##    ggplot2          * 3.3.6    2022-05-03 [2] CRAN (R 4.2.0)
##    ggraph             2.0.5    2021-02-23 [2] CRAN (R 4.2.0)
##    ggrepel            0.9.1    2021-01-15 [2] CRAN (R 4.2.0)
##    GlobalOptions      0.1.2    2020-06-10 [2] CRAN (R 4.2.0)
##    globals            0.15.0   2022-05-09 [2] CRAN (R 4.2.0)
##    glue               1.6.2    2022-02-24 [2] CRAN (R 4.2.0)
##    graphlayouts       0.8.0    2022-01-03 [2] CRAN (R 4.2.0)
##    gridExtra          2.3      2017-09-09 [2] CRAN (R 4.2.0)
##    gtable             0.3.0    2019-03-25 [2] CRAN (R 4.2.0)
##    here               1.0.1    2020-12-13 [2] CRAN (R 4.2.0)
##    highr              0.9      2021-04-16 [2] CRAN (R 4.2.0)
##    htmltools          0.5.2    2021-08-25 [2] CRAN (R 4.2.0)
##    httpuv             1.6.5    2022-01-05 [2] CRAN (R 4.2.0)
##    igraph             1.3.1    2022-04-20 [2] CRAN (R 4.2.0)
##    IRanges            2.30.0   2022-04-26 [2] Bioconductor
##    iterators          1.0.14   2022-02-05 [2] CRAN (R 4.2.0)
##    jquerylib          0.1.4    2021-04-26 [2] CRAN (R 4.2.0)
##    jsonlite           1.8.0    2022-02-22 [2] CRAN (R 4.2.0)
##    knitr              1.39     2022-04-26 [2] CRAN (R 4.2.0)
##    labeling           0.4.2    2020-10-20 [2] CRAN (R 4.2.0)
##    later              1.3.0    2021-08-18 [2] CRAN (R 4.2.0)
##    lattice            0.20-45  2021-09-22 [2] CRAN (R 4.2.0)
##    lifecycle          1.0.1    2021-09-24 [2] CRAN (R 4.2.0)
##    listenv            0.8.0    2019-12-05 [2] CRAN (R 4.2.0)
##    magick             2.7.3    2021-08-18 [2] CRAN (R 4.2.0)
##    magrittr           2.0.3    2022-03-30 [2] CRAN (R 4.2.0)
##    MASS               7.3-57   2022-04-22 [2] CRAN (R 4.2.0)
##    Matrix             1.4-1    2022-03-23 [2] CRAN (R 4.2.0)
##    matrixStats        0.62.0   2022-04-19 [2] CRAN (R 4.2.0)
##    memoise            2.0.1    2021-11-26 [2] CRAN (R 4.2.0)
##    mgcv               1.8-40   2022-03-29 [2] CRAN (R 4.2.0)
##    microbiome         1.18.0   2022-04-26 [2] Bioconductor
##    microViz         * 0.9.1    2022-05-16 [1] Github (david-barnett/microViz@07f1fdb)
##    mime               0.12     2021-09-28 [2] CRAN (R 4.2.0)
##    multtest           2.52.0   2022-04-26 [2] Bioconductor
##    munsell            0.5.0    2018-06-12 [2] CRAN (R 4.2.0)
##    nlme               3.1-157  2022-03-25 [2] CRAN (R 4.2.0)
##    parallelly         1.31.1   2022-04-22 [2] CRAN (R 4.2.0)
##    patchwork          1.1.1    2020-12-17 [2] CRAN (R 4.2.0)
##    permute            0.9-7    2022-01-27 [2] CRAN (R 4.2.0)
##    phyloseq         * 1.40.0   2022-04-26 [2] Bioconductor
##    pillar             1.7.0    2022-02-01 [2] CRAN (R 4.2.0)
##    pkgbuild           1.3.1    2021-12-20 [2] CRAN (R 4.2.0)
##    pkgconfig          2.0.3    2019-09-22 [2] CRAN (R 4.2.0)
##    pkgload            1.2.4    2021-11-30 [2] CRAN (R 4.2.0)
##    plyr               1.8.7    2022-03-24 [2] CRAN (R 4.2.0)
##    png                0.1-7    2013-12-03 [2] CRAN (R 4.2.0)
##    polyclip           1.10-0   2019-03-14 [2] CRAN (R 4.2.0)
##    prettyunits        1.1.1    2020-01-24 [2] CRAN (R 4.2.0)
##    processx           3.5.3    2022-03-25 [2] CRAN (R 4.2.0)
##    promises           1.2.0.1  2021-02-11 [2] CRAN (R 4.2.0)
##    ps                 1.7.0    2022-04-23 [2] CRAN (R 4.2.0)
##    purrr            * 0.3.4    2020-04-17 [2] CRAN (R 4.2.0)
##    R6                 2.5.1    2021-08-19 [2] CRAN (R 4.2.0)
##    RColorBrewer       1.1-3    2022-04-03 [2] CRAN (R 4.2.0)
##    Rcpp               1.0.8.3  2022-03-17 [2] CRAN (R 4.2.0)
##    RCurl              1.98-1.6 2022-02-08 [2] CRAN (R 4.2.0)
##    registry           0.5-1    2019-03-05 [2] CRAN (R 4.2.0)
##  P remotes            2.4.2    2021-11-30 [?] CRAN (R 4.2.0)
##    reshape2           1.4.4    2020-04-09 [2] CRAN (R 4.2.0)
##    rhdf5              2.40.0   2022-04-26 [2] Bioconductor
##    rhdf5filters       1.8.0    2022-04-26 [2] Bioconductor
##    Rhdf5lib           1.18.0   2022-04-26 [2] Bioconductor
##    rjson              0.2.21   2022-01-09 [2] CRAN (R 4.2.0)
##    rlang              1.0.2    2022-03-04 [2] CRAN (R 4.2.0)
##    rmarkdown          2.14     2022-04-25 [2] CRAN (R 4.2.0)
##    rprojroot          2.0.3    2022-04-02 [2] CRAN (R 4.2.0)
##    rstudioapi         0.13     2020-11-12 [2] CRAN (R 4.2.0)
##    Rtsne              0.16     2022-04-17 [2] CRAN (R 4.2.0)
##    S4Vectors          0.34.0   2022-04-26 [2] Bioconductor
##    sass               0.4.1    2022-03-23 [2] CRAN (R 4.2.0)
##    scales             1.2.0    2022-04-13 [2] CRAN (R 4.2.0)
##    seriation        * 1.3.5    2022-03-28 [2] CRAN (R 4.2.0)
##    sessioninfo        1.2.2    2021-12-06 [2] CRAN (R 4.2.0)
##    shape              1.4.6    2021-05-19 [2] CRAN (R 4.2.0)
##    shiny            * 1.7.1    2021-10-02 [2] CRAN (R 4.2.0)
##    stringi            1.7.6    2021-11-29 [2] CRAN (R 4.2.0)
##    stringr            1.4.0    2019-02-10 [2] CRAN (R 4.2.0)
##    survival           3.3-1    2022-03-03 [2] CRAN (R 4.2.0)
##    testthat           3.1.4    2022-04-26 [2] CRAN (R 4.2.0)
##    tibble             3.1.7    2022-05-03 [2] CRAN (R 4.2.0)
##    tidygraph          1.2.1    2022-04-05 [2] CRAN (R 4.2.0)
##    tidyr              1.2.0    2022-02-01 [2] CRAN (R 4.2.0)
##    tidyselect         1.1.2    2022-02-21 [2] CRAN (R 4.2.0)
##    TSP                1.2-0    2022-02-21 [2] CRAN (R 4.2.0)
##    tweenr             1.0.2    2021-03-23 [2] CRAN (R 4.2.0)
##    usethis            2.1.5    2021-12-09 [2] CRAN (R 4.2.0)
##    utf8               1.2.2    2021-07-24 [2] CRAN (R 4.2.0)
##    vctrs              0.4.1    2022-04-13 [2] CRAN (R 4.2.0)
##    vegan              2.6-2    2022-04-17 [2] CRAN (R 4.2.0)
##    viridis            0.6.2    2021-10-13 [2] CRAN (R 4.2.0)
##    viridisLite        0.4.0    2021-04-13 [2] CRAN (R 4.2.0)
##    withr              2.5.0    2022-03-03 [2] CRAN (R 4.2.0)
##    xfun               0.31     2022-05-10 [2] CRAN (R 4.2.0)
##    xtable             1.8-4    2019-04-21 [2] CRAN (R 4.2.0)
##    XVector            0.36.0   2022-04-26 [2] Bioconductor
##    yaml               2.3.5    2022-02-21 [2] CRAN (R 4.2.0)
##    zlibbioc           1.42.0   2022-04-26 [2] Bioconductor
## 
##  [1] /Users/david/Documents/my_projects/evomics/evomics-material/renv/library/R-4.2/aarch64-apple-darwin20
##  [2] /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library
## 
##  P ── Loaded and on-disk path mismatch.
## 
## ──────────────────────────────────────────────────────────────────────────────